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1 Introduction 

A major challenge for molecular modeling consists in optimizing the unlike interac- 
tion potentials. A broad study on fluid mixtures [1] recently showed that among the 
variety of combination rules that were proposed in the past, none is clearly superior. 
In many cases, all are suboptimal when accurate predictions of properties like the 
mixture vapor pressure are needed. The well known Lorentz-Berthelot rule performs 
quite well and can be used as a starting point. If more accurate results are required, it 
is often advisable to adjust the dispersive interaction energy parameter which leads 
to very favorable results [1, 2, 3, 4, 5]. 

A similar approach should be followed for effective pair potentials acting be- 
tween fluid particles and the atoms of a solid wall. They can only be reliable if 
fluid-wall contact effects are taken into account, e.g., by adjusting unlike parameters 
to contact angle measurements [6]. Teletzke et al. [7] used continuum methods to 
examine the dependence of wetting and dewetting transitions on characteristic size 
and energy parameters of the fluid-wall dispersive interaction. MD simulation can be 
applied for the same purpose, leading to a consistent molecular approach. 

On the molecular level, the precise position of the vapor-liquid phase boundary 
is defined by a cluster criterion. Many different criteria are known and it is not im- 
mediately obvious which of them leads to the most accurate results [8]. In nanoscopic 
systems, minute absolute differences can lead to comparably large relative devia- 
tions. Therefore, the viability of several criteria is compared in the present study 
with the purpose of excluding errors due to an inaccurate detection of the interface. 



Author to whom correspondence should be addressed: Prof. Dr.-Ing. habil. J. Vrabec. 
E-mail: jadran.vrabec@upb.de. 



2 Horsch et al. 



If the cohesion of the liquid phase is partly due to hydrogen bonds, successful 
molecular models for pure fluids can often be developed on the basis of an ab ini- 
tio study of the charge distribution as well as the equihbrium position of the nuclei. 
This sterically reaUstic approach, combined with adjusting the Lennard-Jones (LJ) 
potential parameters to vapor-liquid equilibrium ( VLE) data, leads to empirical mod- 
els that correctly reproduce and predict thermophysical fluid properties over a wide 
range of conditions [9]. The present work applies this approach to mixtures contain- 
ing hydrogen bonding components. Often potential parameters determined for one 
fluid carry over to a derivative with different substituents, opening the possibihty of 
creating generic molecular models. Such a model is presented for benzyl alcohol. 

The following pubUcations in peer-reviewed international journals contribute to 
the present project: 

• Schnabel, T, Vrabec, J. & Hasse, H. Molecular simulation study of hydrogen 
bonding mixtures and new molecular models for mono- and dimethylamine. 
Fluid Phase Equilib. 263: 144-159(2008). 

• Eckl, B., Vrabec, J. & Hasse, H. An optimized molecular model for ammonia. 
Mol. Phys. 106: 1039-1046(2008). 

• Eckl, B., Vrabec, J. & Hasse, H. Set of molecular models based on quantum 
mechanical ab initio calculations and thermodynamic data. J. Phys. Chem. B 112: 
12710-12721 (2008). 

• Vrabec, J., Huang, Y.-L. & Hasse, H. Molecular models for 267 binary mixtures 
validated by vapor-liquid equihbria: a systematic approach. Fluid Phase Equilib. 
279: 120-135 (2009). 

• Huang, Y.-L., Vrabec, J. & Hasse, H. Prediction of ternary vapor-liquid equilibria 
for 33 systems by molecular simulation. Fluid Phase Equilib. 287: 62-69 (2009). 

• Horsch, M., Heitzig, M., Dan, C, Harting, J., Hasse, H., Fischer, J. & Vrabec, J. 
Contact angle dependence on the fluid wall dispersive energy. In preparation. 

It would exceed the scope of the present report to give a full exposition of these 
articles. Instead, a few points are emphasized and arranged as follows: Firstly, mix- 
ture properties are explored for binary systems containing hydrogen bonding com- 
ponents. Secondly, vapor-liquid interface cluster criteria and contact angles are dis- 
cussed and remarks on computational details are given. Finally, a sterically accurate 
generic model for benzyl alcohol is introduced and evaluated. 

2 Fluid mixtures with hydrogen bonding components 

Vapor- liquid equihbria of 3 1 binary mixtures consisting of one hydrogen bonding 
and one non-hydrogen bonding component were studied. All models are of the rigid 
united-atom multi-center LJ type with superimposed electrostatic sites in which hy- 
drogen bonding is described by partial charges. The hydrogen bonding components 
of the studied binary mixtures are: monomethylamine (MMA) and dimethylamine 
(DMA), methanol, ethanol and formic acid. The non-hydrogen bonding components 
are: neon, argon, krypton, xenon, methane, oxygen, nitrogen, carbon dioxide, ethyne. 
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ethene, ethane, propylene, carbon monoxide, diflourodichloromethane (R12), tetra- 
flouromethane (R14), diflourochloromethane (R22), diflouromethane (R32), 1,1,1,2- 
tetraflouroethane (R134a) and 1,1-diflouroethane (R152a). 

To obtain a quantitative description of the mixture vapor-liquid equilibria, one 
state independent binary interaction parameter was adjusted to a single experimental 
data point of either the vapor pressure or the Henry's law constant. Throughout, 
excellent predictions were found at other state points, i.e., at other compositions 
or temperatures as well as for the Henry's law constant, if it was adjusted to the 
vapor pressure, or vice versa. Figures 1 and 2 show methane + methanol as a typical 
example. 




0.00 0.05 0.10 0.96 1.00 

Xr^^A I mol mol"^ 



Fig. 1. Simulation data and vapor-liquid equilibria of methanol + methane: simulation data 
(•) and experimental data (o) at 310 K [10]; simulation data (A) and experimental data (A) at 
338.15 K [11]. 

Furthermore, a set of 78 pure substances from prior work [22, 23] was taken 
to systematically describe all 267 binary mixtures of these components for which 
relevant experimental VLB data are available. Again, per binary system, the single 
state independent binary interaction parameter in the energy term was adjusted to 
only one experimental value of the vapor pressure. The unlike energy parameter 
was thereby altered usually by less than 5 % from the Berthelot rule. The mixture 
models were validated regarding the vapor pressure at other state points as well as 
the dew point composition, which is a fully predictive property in this context. In 
almost all, i.e., 97 % of the cases, the molecular models give excellent predictions of 
the mixture properties. Compared to other works in the Uterature, this is by far the 
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Fig. 2. Henry's law constant of methane in methanol: simulation data (•) and experimental 
data [12, 13, 14, 15, 16, 17, 18, 19, 20, 21] (o). 



largest investigation in this direction. It was facilitated by the extensive computing 
equipment at the High Performance Computing Center Stuttgart. 

In the next step, all 33 ternary mixtures of these 78 components for which exper- 
imental VLE data is available were studied by molecular simulation. No adjustment 
to ternary data was carried out at all so that the calculations were strictly predictive. 
By comparing to experimental data, it was found that these predictions are very reh- 
able as there was practically always an excellent match. As an example, Fig. 3 shows 
the ternary system consisting of methane, ethane, and carbon dioxide. Again, the 
computational effort was substantial, publications in the literature by other groups in 
this field typically cover one to two mixtures only. 



3 Vapor-liquid interface cluster criteria 

A suitable cluster criterion should achieve two goals: on the one hand, it needs to 
distinguish the bulk Uquid and the bulk vapor successfully in every case - even when 
the vapor is supersaturated, such as in the vicinity of a droplet, or the liquid is under- 
saturated, such as in the vicinity of a bubble. On the other hand, the cluster criterion 
should also minimize noise fluctuations of the detected clusters to emphasize the 
signal. 

The following criteria for carbon dioxide were compared for this purpose using 
a rigid two-center LJ plus point quadrupole (2CLJQ) model [22]: 
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Mole fraction C^Hg 



Fig. 3. Ternary vapor-liquid equilibrium phase diagram for the mixture CH4 + CO2 + CaHg 
at 230 K and 4650 kPa: simulation data (•), experimental data (+) of Knapp et al. [24], and 
Peng-Robinson equation of state ( — ). 

• Stillinger [25]: all molecules with a distance of rst or less from each other are 
liquid and belong to the same cluster (i.e., the same liquid phase or the same 
droplet). The Stillinger radius was set to rst = 5.7 A for the present simulations. 

• Ten Wolde and Frenkel [26]: all molecules with at least k neighbors within a 
distance of rst belong to the liquid. They belong to the same cluster as all other 
liquid molecules within a distance of rst. 

• Arithmetic mean, k neighbors: a molecule is liquid if the density in the sphere 
containing its k nearest neighbors exceeds (p' + p")/2, where p' amd p" are 
the saturated liquid and vapor density, respectively. The molecule belongs to the 
same cluster as aU other liquid molecules within the radius ri^, which defines a 
sphere with the volume occupied by ^ + 1 molecules at a density of (p' + p")/2. 

• Geometric mean, k neighbors: analogous, with a density threshold of (p'p")^/^. 

The simulations were conducted with the massively parallel MD program Isl mar- 
dyn (the precursor implementation Isl moldy is described by Bernreuther and Vrabec 
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[27]). Figure 4 shows that all of the discussed criteria are applicable. The Stillinger 
criterion and the geometric mean density criterion with two neighbors lead to the 
best results. It should be noted that at high temperatures, i.e., near the critical point, 
the Stillinger criterion becomes less reliable in distinguishing the liquid from a su- 
persaturated vapor than the geometric mean density criterion. 



4 Contact angle dependence on the fluid-wall interaction 

In cases where no hydrogen bonds are formed between the wall and the fluid, vapor- 
solid and liquid-solid interfacial properties mainly depend on the fluid-wall disper- 
sive interaction, even for hydrogen bonding fluids. The truncated and shifted LJ 
(TSLJ) potential [28] 




„-6^ 



(1) 



with a cutoff radius of rc = 2.5 a accurately reproduces the dispersive interaction if 
adequate values for the size and energy parameters a and e are specified. Due to the 
relatively small cutoff radius, all computations are accelerated, while the descriptive 




droplet size in number of molecules 



Fig. 4. Average number of molecules evaporating during a detection interval of 30 fs from 
droplets in a supersaturated vapor of carbon dioxide with T = 237 K, p = 1.98 mol/1, and A' 
= 1,050,000, determined according to various cluster criteria: Stillinger (o), ten Wolde and 
Frenkel with k = 4 (•), arithmetic mean with two (A) and eight neighbors (V) as well as 
geometric mean with two neighbors (□). 
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power of the full LJ potential (without a cutoff) is retained even for systems with 
phase boundaries [29]. 

For the present series of contact angle simulations, the TSLJ model - with the 
potential parameters for methane, a = 3.7241 A and e/k = 175.06 K - was studied, 
extending a previous investigation of the vapor-liquid interface of the TSLJ fluid 
[29]. The TSLJ potential with the size and energy parameters o' = O as well as 

e' = ^e, (2) 

was also applied for the unlike interaction between the fluid molecules and the wall 
atoms, with the same cutoff radius as for the fluid. The wall was modeled as a system 
of coupled harmonic oscillators with different spring constants for transverse and 
longitudinal motion, adjusted to simulation results for graphite that were obtained 
with a rescaled version of the Tersoff [30] potential. 

The simulations were carried out with the Isl mardyn program. Vapor and liquid 
were independently equilibrated in homogeneous simulations for 10 ps. This was 
followed by 200 ps of equilibration for the combined system, i.e., a liquid meniscus 
surrounded by vapor, with a graphite wall consisting of four to seven layers, cf. 
Fig. 5. A periodic boundary condition was applied to the system, leaving a channel 
with a diameter of 27 a between the wall and its periodic image. The contact angle 
was determined from the density profiles by averaging over at least 800 ps after 
equilibration. 




Fig. 5. Simulation snapshots for a reduced fluid-wall dispersive energy ^ of 0.07 (left) and 

0. 13 (right) at a temperature of 0.88 t,/k. 

Liquid menisci between graphite walls were simulated for a reduced fluid-wall 
dispersive energy ^ between 0.07 and 0.16 at temperatures of 0.73, 0.88, and 1 z/k. 
Note that the triple point temperature of the TSLJ fluid is about 0.65 t/k [31] while 
the critical temperature is = 1.0779 t/k [29], so that the entire regime of stable 
vapor-liquid coexistence was covered. 

High values of ^ correspond to a strong attraction between fluid and wall com- 
poents, leading to a contact angle smaller than 90°, i.e., to partial (■& > 0°) or full 
("& = 0°) wetting of the surface. For a higher fluid- wall dispersive energy, the extent 
of wetting increases, cf. Fig. 5. The transition from obtuse to acute contact angles, 

1. e., i3-(r,^o) 90°, occurs at a temperature independent value ~ 0.113 of the 
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fluid-wall dispersive energy, as can be seen in Fig. 6. Furthermore, the symmetry law 

cosiJ(r,^o -A^) - -cosiJ(r,^o + A^), (3) 

is valid over the whole relevant range of temperatures and magnitudes of the fluid- 
wall dispersive energy. Figure 6 also shows that there is a narrow range of t, values 
that lead to the formation of a contact angle, as opposed to total dewetting or wetting. 
As the temperature decreases and the vapor-hquid surface tension "/{T) increases, the 
contact angle approaches 90°. 




0.050 0.075 0.100 0.125 0.150 0.175 
fluid-wall interaction energy parameter ^ 

Fig. 6. Simulation results and correlation for the contact angle in dependence of the reduced 
fluid-wall dispersive energy t, at temperatures of 0.73 (diamonds and solid line), 0.88 (circles 
and dotted line) as well as 1 (squares and dashed line) z/k. 



5 A sterically accurate generic benzyl alcohol model 

Benzyl alcohol (C6H5-CH2OH) is widely used as a solvent for paints and inks. How- 
ever, it is classified as a harmful substance (Xn) and should not be inhaled, nor used 
at high temperatures where it exhibits a high vapor pressure. 

The basis for a new rigid molecular model of benzyl alcohol was determined 
by ab initio calculations with the GAMESS (US) quantum chemistry package [32], 
obtaining the quadrupole moment as well as the equilibrium positions of the nuclei 
as illustrated in Fig. 7. 

Further potential parameters were taken from accurate empirical molecular mod- 
els for related fluids, leading to a sterically accurate generic model, cf. Tab. 1, that 
can be used as a starting point for parameter optimization. In particular, point charges 
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as well as a and e parameters for the hydroxyl group were taken from the ethanol 
model of Schnabel et al. 133]. Moreover, the o and e values of the corresponding LJ 
interaction site of the Merker et al. [34] cyclohexanol model were used for the CH2 
group, while the LJ parameters for the CH and C centers were set according to the 
Huang et al. [35] models of benzene and phosgene, respectively. 



Table 1. Coordinates and parameters of the LJ sites, and the point charges, and the point 
quadrupole for the present benzyl alcohol model. Bold characters indicate represented atoms. 



Interaction 
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Z/kB 
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A 
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K 
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CHP 


0.695 


-0.037 


2.218 


3.247 


88.97 






CH'-™ 


-0.753 
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3.247 
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0.189 
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3.247 
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-0.285 
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-0.497 
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CH2 


-1.160 
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+0.2556 




OH 
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-0.6971 
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Fig. 7. Coordinates of the LJ sites for the present benzyl alcohol model. 
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Fig. 8. Vapor pressure of benzyl alcohol according to a correlation [36] based on experiments 
( — ) as well as present molecular simulation data (o). 

VLE properties of the generic model were calculated using the ms2 program, 
leading to an overall satisfactory first approximation, considering that all a posteriori 
adjustments were absent. The vapor pressure is more accurate at high temperatures, 
cf. Fig. 8. 

6 Computing performance 

The scalability of the Isl mardyn program was measured on the cacau supercomputer 
for simulation scenarios involving methane, represented by the TSLJ fluid, as well 
as graphite, modeled by a rescaled version of the Tersoff [30] potential. MPI par- 
allehzation was appUed according to a spatial domain decomposition scheme with 
equally sized cuboid subdomains and a cartesian topology based on linked cells [27]. 

Often the best solution is an isotropic decomposition that minimizes the surface 
to volume ratio of the spatial subdomains. For the simulation of homogeneous sys- 
tems, this approach is quite efficient. That is underlined by the weak and strong 
scaling behavior of Isl mardyn for typical configurations, shown in Fig. 9 (left), in 
cases where supercritical methane ('fluid') at a density of 10 mol/1 and soUd graphite 
('wall') were considered with a system size of up to 4,800,000 interaction sites, rep- 
resenting the same number of carbon atoms and methane molecules here. Graphite 
simulations, containing only carbon atoms, scale particularly well, due to a favor- 
able relation of the delay produced by communication between processes to the con- 
current parts, i.e., the actual intermolecular interaction computation, which is much 
more expensive for the Tersoff potential than the TSLJ potential. 
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1 10 100 50 100 150 

number of processes number of processes 

Fig. 9. Left: Total CPU time, i.e., execution time multiplied with the number of parallel pro- 
cesses, per time step and interaction site for weak scaling with 3,000 (dashed lines / circles) 
and 32,000 (dotted lines / triangles) interaction sites per process as well as strong scaling with 
450,000 interaction sites (solid lines / bullets), using isotropic spatial domain decomposition. 
Right: speedup for a system of liquid methane between graphite walls with 650,000 inter- 
action sites, where isotropic (circles) and channel geometry based (triangles) spatial domain 
decomposition was used; the solid line represents linear speedup. 

The simulation of combined systems, containing both fluid and solid interaction 
sites, is better handled by a channel geometry based decomposition scheme, where 
an approximately equal portion of the wall and a part of the fluid is assigned to each 
process, cf. Fig. 9 (right). In the general case, where spatial non-uniformities do not 
match any cartesian grid, a flexible topology has to be used. An approach based on 
A:-dimensional trees [37], implemented in a version of Isl mardyn, showed clearly 
improved results with respect to the scaling of inhomogeneous systems. 

7 Conclusion 

The intermolecular interaction of small hydrogen bonding molecules like mono- and 
dimethylamine can be described by simple LJ based united-atom molecular models 
with point charges. Such computationally efficient models were applied to binary 
mixtures with non-hydrogen bonding components regarding VLB properties. Ac- 
curate predictions, covering a broad range of temperatures and compositions, were 
obtained, regardless whether the state independent binary interaction parameter was 
adjusted to Henry's law constant or vapor pressure. 

For a two-center LJ plus quadrupole model of carbon dioxide, a comparison of 
cluster criteria with the purpose of accurately detecting the vapor-liquid phase bound- 
ary gave overall support to the geometric mean density criterion applied to the sphere 
consisting of a molecule and its two nearest neighbours. The Stillinger criterion was 
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found to be particularly adequate at low temperatures. For the TSLJ fluid, the con- 
tact angle formed between the vapor-liquid interface and a wall was determined by 
canonical ensemble MD simulation while the mangitude of the dispersive fluid-wall 
interaction was varied. Over the whole temperature range between triple point and 
critical point, the contact angle dependence on the fluid- wall dispersive energy obeys 
a simple symmetry law. 

A sterically realistic model for benzyl alcohol was presented, showing within 
the framework of LJ sites with point polarities and electric point charges that the 
generic molecular modeling approach can lead to a good starting point for parameter 
optimization with respect to VLE properties such as the vapor pressure. 

The scalability of the Isl mardyn program was assessed and found to be ac- 
ceptable. MD simulations of methane confined between graphite walls with up to 
4,800,000 interaction sites, i.e., carbon atoms and methane molecules, were con- 
ducted to demonstrate the viability of the program. 

The authors would hke to thank Martin Bernreuther, Martin Buchholz, Domenic 
Jenz, and Christoph Niethammer for contributing to the Isl mardyn code, Stephan 
Deublein, Bernhard Eckl, Gimmy Fernandez Ramirez, Gabriela Guevara Carrion, 
and Sergei Lishchuk for contributing to the ms2 code, loannis Pitropakis for per- 
forming simulation runs of binary mixtures, as well as CaUn Dan, Franz Gahler, Se- 
bastian Grottel, Jens Harting, Martin Hecht, Ralf Kible, and Guido Reina for frank 
discussions and their persistent support. 
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